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1.  INTRODUCTION 


It  is  generally  accepted  that  the  transition  of  laminar  flow  to  turbu¬ 
lent  flow  in  channels  and  boundary  layers  is  preceded  by  amplification 
of  smal 1 -ampl i tude  disturbances  super-imposed  on  the  steady  laminar 
flow.  In  the  early  stages,  nonlinear  effects  are  negligible,  and  fur¬ 
thermore,  if  the  slow  variation  of  undisturbed  flow  with  distance  along 
the  flow  is  neglected  -  quasi -one-dimensional  flow  -  the  small  disturb¬ 
ances  may  be  decomposed  into  Fourier  components  In  time  and  streamwise 
distance.  The  Fourier  components  satisfy  the  homogeneous  Orr-Sommerfeld 
equation,  whose  eigenvalues  determine  the  stability  of  the  Fourier  com¬ 
ponents. 

The  results  from  linear  stability  theory  are  applied  in  the  prediction 
of  transition  by  Smith  (Smith  and  Gamberonl,  1956;  and  Jaffe,  Okamura 
and  Smith,  1970)  and  by  Mack  (1975).  The  Smith  method  is  based  on  the 
correlation  of  the  amplitude  ratio  computed  by  linear  theory  with 
experimental  data  on  the  transition  Reynolds  number.  Transition  Is 
assumed  to  take  place  when  the  computed  amplitude  for  the  most  amplified 
wave  reaches  eN  with  N  ranging  from  8  to  11.  The  Mack  method,  on  the 
other  hand,  is  based  on  the  magnitude  of  the  computed  amplitude.  The 
initial  amplitude  of  a  disturbance  at  the  critical  Reynolds  number  Is 
estimated  from  the  power  spectrum  of  free-stream  turbulence,  and  the 
subsequent  amplitude  Is  computed  based  on  linear  stability  theory.  Tran¬ 
sition  Is  assumed  to  take  place  when  the  computed  amplitude  reaches  a 
certain  reference  value.  Mack  applied  the  method  to  predict  the  effect 
of  freestream  turbulence  level  on  the  transition  Reynolds  number  and 
obtained  good  agreement  with  the  measured  data. 

These  prediction  methods  apply  linear  stability  theory  which  one  fully 
recognizes  does  not  apply  when  the  disturbance  amplitudes  are  no  longer 
very  small,  and  one  would  hope  that,  if  nonlinear  effects  are  Included 
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1n  the  computation,  one  could  predict  transition  with  better  precision 
and  confidence.  Inclusion  of  nonlinear  effects,  however,  magnifies  the 
complexity  by  a  large  factor.  We  do  not  have  the  exact  analytical  solu¬ 
tion  and  are  forced  to  employ  some  approximation  methods.  Because  the 
problem  is  nonlinear,  we  cannot  superimpose  elementary  solutions  to 
generate  a  new  solution.  Every  plausible  mechanism  has  to  be  studied 
individually.  Though,  at  present,  the  consensus  among  the  people  who 
are  investigating  boundary-1 ayer  transition  Is  that  the  three-dimen¬ 
sional  effects  are  an  essential  element  in  the  nonlinear  regime  leading 
to  transition.  The  nonlinear  stability  of  two-dimensional  disturbances 
has  been  Investigated  by  many  people,  beginning  with  the  pioneering  work 
of  Stuart  and  Watson  (I960),  which  is  based  on  the  expansion  in  small 
amplitude  and  the  first  approximation  is  the  linear  stability  theory. 
The  nonlinear  stability  theory  Is  In  reasonably  good  shape  for  confined 
flows  that  are  weakly  unstable  in  the  linear  regime,  for  example,  two- 
dimensional  Poiseullle  flow  (Itoh,  1977),  The  theory,  however,  has  not 
been  worked  out  for  boundary-layer  flows,  because  the  treatment  of  the 
effect  on  the  mean  flow  Is  more  subtle  in  boundary-layer  flows  than  in 
the  flow  between  parallel  walls. 

Quite  apart  from  this  analytical  approach,  numerical  integration  of  the 
Navier-Stokes  equation  was  undertaken  to  study  the  stability  of  large 
amplitude  disturbances  (Fasel,  1976;  Murdock,  1977).  Between  the  two 
computations,  however,  there  was  some  disagreement  on  the  relative  phase 
between  the  fundamental -frequency  component  and  the  second-harmonic  com¬ 
ponent. 

The  present  investigation  was  undertaken  to  establish  a  nonlinear  sta¬ 
bility  theory  for  two-dimensional  disturbances  In  boundary-layer  flows. 
The  results  of  the  present  study  serve  as  a  check  on  the  numerical 
Integration  of  the  Navier-Stokes  equation.  This  type  of  approach  Is 
expected  to  give  a  good  representation  of  the  early  development  of  the 
nonlinear  effect  for  a  growing  disturbance  and,  presumably,  can  provide 
further  insight  into  the  boundary-layer  transition  process. 
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The  basic  assumptions  made  In  the  formulation  are: 

(1)  The  amplitude  of  disturbances  Is  small. 

(2)  The  rate  of  change  of  the  undisturbed  boundary  layer  with  stream- 
wise  distance  Is  small. 

(3)  The  rate  of  amplification  of  disturbance  is  small. 

The  first  assumption  allows  us  to  expand  the  solution  In  powers  of  small 
amplitude.  In  the  present  work  we  retained  terms  to  the  third  power  of 
the  amplitude,  but  we  did  not  calculate  the  third  harmonic  in  our  analy¬ 
sis  since  it  does  not  affect  the  form  of  the  evolution  equation  to  the 
order  we  considered.  The  third  assumption  together  with  the  second 
assumption  assures  that  the  mean  flow  (disturbed  as  well  as  undisturbed) 
variation  with  the  streamwlse  distance  Is  small,  so  that  the  disturbed 
mean  flow  remains  boundary-layer  like. 


Mathematical  formulation  Is  presented  in  Section  2,  the  detailed  discus¬ 
sion  of  numerical  method  in  Section  3,  and  the  discussion  of  results  In 
Section  4. 
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2.  FORMULATION  OF  NONLINEAR  STABILITY  PROBLEM 

In  this  section  we  formulate  a  nonlinear  stability  theory  for  two-dimen¬ 
sional  disturbances  in  an  incompressible  boundary  layer  along  a  flat 
plate.  Including  the  effect  of  non-parallel  ness  of  the  mean  flow.  Our 
starting  point  Is  the  Navler-Stokes  equation  written  In  terms  of  the 
streamfunction 

(V2*)t  +  yv2<Mx  -  *x(v2*)y  =  (2.1) 

I 

with  boundary  conditions 

/ 

♦  s  ♦ y  =  0  at  y  =  0,  *  +  1  as  y  ♦ 

Here,  x,y,t  and  <|>  are  made  dimensionless  by  dividing  distances  by  a  re- 

V 

fere^ce  length  i  (comparable  with  the  boundary-layer  thickness),  time  by 
fc/Ua,  (UM  freestream  velocity)  and  streamfunction  by  Ujt.  R  denotes  the 
Reynolds  number  U J-/v. 

The  undisturbed  flow  changes  very  slowly  with  x,  and  we  denote  the  rate 
of  change  by  a  small  parameter  u.  We  introduce  a  slow  variable 

X  =  ux  (2.2) 

and  then  the  basic  flow  is  given  by 

?  *  ?(X,y) 

so  that 

Ty  =0(1)  ,  =  U?X  =  0(w) 


(2.3) 


Furthermore,  we  assume  that  the  perturbation  depends  on  x  and  t  (fast 
variables)  only  through  the  complex  phase  angle  0,  defined  by 


0x  =  k  *  et  =  *“  (2.4) 

The  wave  number  k,  however,  depends  on  the  slow  variable  X.  We  assume  w 
to  be  real,  but  k  complex  with  a  small  Imaginary  part;  namely  we  have  a 
spatial  stability  problem  with  a  weak  amplification  rate. 

Then,  the  streamf unction  assumes  the  following  form: 

*  ■  *(0)(X,y)  +  Xj  j*(n)(X,y)  e1n  +  ?{n)(X,y)  e  ^  j  (2.5) 

where  C)  denotes  the  complex  conjugate.  Substituting  the  above  expres¬ 
sion  into  equation  (2.1)  and  collecting  the  terms  with  equal  values  of 
n,  we  obtain  the  following  infinite  set  of  equations: 

MEAN  FLOW: 


irA.2  *(0)-  <j4°)+  40)  A  f(0) 


Wo 


y  o  yX  rX  o  y 


00 

=  i-Xje  2  1  *(-ink+u3x)  An^n)  +  ^{ink+uS^A^"* 


(ink+uav)  •  An  -  (-1nk+u3v)  •  a, 


n  y 


,  M 

n  y  j 


(2.6) 


t 
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FUNDAMENTAL: 

I  A2  iajA  <j/°'(1k+y3y)  A  i/1* 

K  1  1  y  X  i 

+  ( )  (ik+p3x)  *(1)-  u(aq^0))«^1)+  w^0)  a^1* 

00 

=Xje  1  |^1)(-fmit+M3x)Am^n,)+^m)[l(nH-l)k+w3x]Am+1*(nB+1) 

-  [i(n,+l)k+u3x]*{m+1).Am^m)  -  (-1me+li3x)?(ra).Am+1^m+1)|  (2.7) 

HIGHER  HARMONICS:  n  =  2 

^  A24»(nUina)An^(n)-^0)(1nk+M3x)An^n)+(Ao^(0))(1nk+w3x)ip(n) 

•  m( Ao'pX°>)’iyn>+UTf'x0>  Vy"* 

00 

=  Zj  |^n'm^1mk+M3x)Am^mL[i(n-m)k+w3x]^{n'm).Am^m)| 

°°  2  0 

♦Ee  1  }*in™,(-i«*t»\)4„»<")+*'*,[1(n«>)ktu8Y]A  *(n™> 

1  f  J  Amy  A  n '*in 

-[1(n+m)k+M3x]^n+m).  Vym)-("'lnk+u3x)?(m)*  An+ij)^n+m)  |  (2.8) 


In  these  equations,  the  following  notations  are  used: 


3X  =  partial  differentiation  with  respect  to  the  slow  variable  X. 

3y  =  partial  differentiation  with  respect  to  y. 

An  =  3^  +  (ink+u3x)2 

=  32  -  n2k2  +  uin(2k3x+kx)  +  p23^ 

Mean  Flow  Solution  -  Lowest  Order: 

To  the  lowest  order,  the  mean-flow  equation  Is 

*-»y  Sy  3y*=0 

for  =  J  +...  .  Integrating  once  with  respect  to  y,  we  obtain 

k  *m-\%  *  Wy  '  p,x)  <2-9) 

which  we  recognize  as  the  usual  boundary -layer  equation,  after  identify¬ 
ing  P ( X)  to  be  the  pressure  gradient  Imposed  by  the  flow  external  to  the 

boundary  layer. 

Fundamental  Component  -  Lowest  Order 

The  fundamental  component  equation  to  the  lowest  order  Is  the  Orr- 
Sommerfeld  equation: 

Lu(*(1))  =  £  02-k2)%(1)-  l[(kUy-a,)(32-k2)^(1)-k^y^(1)]=0  (2.10) 
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and  the  solution  is 

*(1)  =  A(X)$(1)(y;X,<o)  (2.11) 

where  A(X){<0(1))  is  the  amplitude  function,  and  satisfies 
=  0.  With  the  boundary  conditions; 

*{1)(0)  «  ♦“>«»  •  *(1)(-)  -  o. 

Thus,  we  obtain  the  eigenvalue  solution 

$(1)  =  *(1)  (y;X,nj)  (2.12a) 

k  =  k(X,«)  (2.12b) 


Second  Harmonic  Component 

The  equation  for  the  second-harmonic  component,  generated  by  the  above 
fundamental  component,  is  then 

l2„(*(2))  s  |  (3j-4k2)%(2)-2i|(kyu))(32-4k2)ip(2)-k?yyy^(2)j 

■  ikA^XjL^V1*  -  +... 

|_  y  yy  yyyj 

Therefore,  to  the  lowest  order,  the  second  harmonic  will  be  of  the  form: 
*(2)  =  A2(X)  <t>( 2) (y ;X,u») 

l2  (»<2>)  -  1kL(1U{1)  -  *{l)M 

ry  yy  y  YyyyJ 


(2.13a) 

(2.13b) 


-9- 


Correctlon  to  Mean  Flow: 

Also,  from  the  lowest-order  fundamental  component,  we  obtain  the  follow¬ 
ing  equation  for  the  correction  to  the  mean  flow, 

u!T  Wlyyyy  vy  VlXyy  vXyy  Vly  VX  V1yyy  vyyy  V1X 


=  i  M  .... 

The  neglected  terms  are  of  the  order  of  u2  and  e4/w. 

Integrating  once  with  respect  to  y,  we  get 

*  ,!»)  .  f  ,!<»  +  j  .  f  ,!«).  f  ,(°l 

HIT  Vyyy  wy  VlXy  VX  Vlyy  vXyVly  vyyVlX 

*iAAe  1  -  k#(1,*J1))  -  4ki^y1>^y1>J+  pi00  +••• 

where  Pj ( X)  again  Is  a  term  due  to  streamwlse  pressure  gradient  and  will 
be  0(m2)  for  a  flat-plate  boundary  layer.  We  assume  |A|2/p>>u2,  and 
let 


•26, 


,i°>  .  | A , 2  e  M°’(X,y) 

Thus,  the  equation  for  becomes 

i  -?  4°^  4°^  -  ?  4°^  +  J  4°*  +  2k. (?  40^-*  4°^) 

R  Jyyy  y  lxy  Yx  Ylyy  rxy  Yly  ry  yix  1'*yyiy  *yyyi  ' 


(2.15a) 


(2.15b) 


In  this  equation,  we  have  neglected  d/d X  *n|A|  compared  with  k^*6jX 
since  the  former  Is  of  the  order  of  u.  Also,  we  have  gone  back  to  the 
original  x  by  multiplying  through  by  m. 
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Corrections  to  Fundamental  Component: 

The  modification  to  the  fundamental  component  consists  of  two  parts: 
one  part  due  to  the  variation  of  and  k  with  X,  and  the  other  part 

due  to  nonlinear  feedback.  By  putting 

,  klA  +  i  \A2A  (2.16) 

and 

x  A*(1)  +  uA*iU  +  | A | 2  a4X)  (2.17) 

we  obtain  the  following  equations  for  the  corrections  and 

L^ul1*)  *  kiMU(1))  +  Q<?,*{1),k)  (2.18) 

Lw(^n)  =  XM<*(1))  +  (2.19) 

where , 

H(^1,)-2uk*(l,+‘yy(^,-3k2#(1))-  k( ♦yy'”k2d>( )  (2.20a) 


Qd'.<>(1),k)^yyx^1)-«(>x(^-k2^1))  -  ?yyx  ^1} 

*  2"k*iU*  '  3k2*i1)>  -  r- k  <♦$  -  k2<tin) 

-  3k?y  t'1'  -  |i  -  3k2*(1))]kx 


(2.20b) 
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and 

+  ^(1)(*yyy  ‘  4k2^2))  -  2k*(2)(*yyy  -  k2?y1}> 

♦  UyJ’  -  k2*(1))  -  k*(JJy  ♦{1)|  (2.20c) 

The  boundary  conditions  for  (n  *  1,2)  are 

♦in(°)  3  ♦JJ,(0)  »  ♦J1)(»)  *  0 

The  equation  LuUj^]*  0  with  above  boundary  conditions  has  nontrivial 
solutions  (elgensolutlons).  Therefore,  If  the  nonhomogeneous  equations 
(2.18)  and  (2.19)  are  to  have  solutions,  the  terms  on  the  right  hand 
side  must  be  orthogonal  to  the  elgensolutlons  of  the  homogeneous  equa¬ 
tion  (hereinafter  this  Is  referred  to  as  the  solvability  condition), 
which  yields: 


[>>* Q  dy 

id - 2_ - 

/  *(1)*M  dy 
o 

/  dy 

A  .  .  1  -2_ - 

/  dy 

o 

111* 

where  ♦  Is  the  solution  of  the  adjoint  problem 


(2.21) 


(2.22) 
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1  (a2-k2)2  ♦(n*.  l[(ki;y-“H32-k2)2  ♦(1)*+  0 


with 


♦  «ll*(0>  *  *  0. 


(2.23) 
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SUMMARY: 


We  summarize  the  mathematical  problems  that  need  to  be  solved: 
(1)  Fundamental -Mode  Eigenvalue  and  Elgensolutlon 

I  U<1>  k  1  =  1  82  k2  2 


(k  ¥  -to)  a*1*  -  k2  -  k  IF  a*1* 

1  o9y  '  fyy  o  9  o9yyy  9 


with  4>^(0)  *  ♦^(O)  *  0  ;  <fr^(y)  +  0  as  y  ♦  • 


Here  Is  the  Orr-Sommerfeld  linear  differential  operator. 


(2)  Adjoint  Orr-Sownerfeld  Problem 


subject  to  the  boundary  conditions,  *  #^*(0)  ■  0, 

♦^*(y)  ♦  0  as  y  ♦  •  and  with  k0  obtained  from  (1). 


/ 
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(3)  Correction  to  Fundamental  Mode  and  Determination  of  Amplification 
Factor  Due  to  Non-Parallel  Effects 


Lu»(*l1)’ko)  *  kl  +  0U(1)) 


with  boundary  conditions  Identical  to^  those  for 
kj  Is  determined  from  the  solvability  condition: 


k(  s  k1(»(U,»ll)’,k<>)  ■  ° 


/V1'*  « 


where 


*  1 


/  dy 

o 

[2u»k  U(l)-3k2*(1))  -  f  4(l)l 

L  o*  y  yy  o  yyy  J 


ko(^}’-ko  ♦,1)) 


and 


QU(l)) 


I ^ 2uik  ♦  uk  *(1)-  f  M 
Y o9x  ox’  vyyy*%  ( 


♦  1?  -  3k2^^^  -  3k  k 

y  |  yyx  JVx  JVox*  j 

1  (ai,  +  7k  aif  3A(1)  ck2k  AD 

'  ff  }4ko*yyX  +  2koX*yy  '  4ko*X  ‘  6kokoX*  , 


♦  fc  .  tT  -  k2^1*) 

*yyx  9y  1vXl,yyy  Vy  ' 


f 
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(4)  Equations  for  and  kQX 


and  kQX  appearing  in  Q  are  obtained  by  solving  the  equation 


■  lHkox*pl 


subject  to  the  boundary  conditions  identical  to  those  for  V 
Of  the  terms  on  the  right  hand  side,  P  Is  defined  by 

JDl 


(1) 


=  ik0  uLV  -  ^ 


xy  Tyy  o 
M  is  as  defined  in  (3). 


xyyy 


(5)  Second  Harmonic  Mode 

L2«'»(21>2kol  ‘  F2(»(1)-k0) 

where  L2w  is  the  Or r- Somme rf el d  operator  with  k0  obtained  from  (1) 

and  F2Ull),k0)  »  -  -  ♦yU4J’]- 

with 

*(2)(0)  =  ^2)(0)  =  0,  *(2)(y)  ♦  0  as  y  +  - 


T 
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(6)  Correction  to  Mean  Flow  Due  to  Fundamental  Mode  and 
Non-Parallel  Effects 

G(*(0))  *  F1(^(1),k0) 


where 


7  ♦<»> 
yx  y 


*  ,<°>  ,  *  ♦«)  .  ?  ,«» 

yyx  y  yx  xyy 


‘  ♦«>  *  2k  ,4  *(01  .  »  ,(0f 

If  9yyy  oi\v,,w  9  v 


yy 


y  y 


') 


F1U<ll,k0)  *  i[k0*<11*y1)-VU)*yU] 


♦  4k0t  UyU| 


subject  to  the  boundary  conditions. 


for  all  x  =  Xj 


For  the  Initial  condition  we  arbitrarily  choose; 
=  0  at  x  =  Xj,  for  all  y. 


^  =  o  at  y  =  0 


JO) 


(as  y  ♦  «) 
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(7)  Correction  of  Fundamental  Mode  and  Determination  of  Landau 
Constant  x,  Due  to  Nonlinear  Effects 


Lj41)»k0) =  xm(*(1))  +  1  N  U(0),*{1),*{2)) 


where  M  and  N  are  defined  below.  The  boundary  conditions  are  iden¬ 
tical  to  those  for 


The  Landau  Constant  X  is  determined  from  the  solvability  condition 
that  the  right-hand- side  of  the  above  equation  has  to  be  orthogonal 
to  Thus,  we  have 

f  N<j>(1)*dy 


f  M*(1)*dy 


where 


+  2k  ?^($^-4k2^2))  -  2k  -  k2$(1)) 

o  y  yy  oT  oy  vvyyy  ovy  ' 


+  k  *(1)(*<2)  -  4k2*(2) 

Ko?  1 9yyy  ^Ko9y 


)  +  M 


oTy  ' ryy 


(1), 


-  k  *{0)M 

oYyyy*  j 

and  M  is  as  defined  in  Problem  (3)  and  kQX  is  determined  by  the 
solvability  condition  to  be 

-fp  *(1)* d y 

U  a  0 _ 

ox  m* 

J  M  r1  dy 
o 


(8)  The  Amplitude  Equation 


A  •  k  A  +  k.A  U  A2  A 
x  o1  1 

subject  to  the  appropriate  Initial  condition  for  A.  The  quantities 
kQ,  kj  and  X  are  determined  in  Problems  (1),  (3)  and  (7),  respec¬ 
tively.  By  expressing  A  in  the  polar  form 

A  =  l  A 1  eU 

and  substituting  into  equation  (2.16),  we  obtain  the  following 
equations  for  the  amplitude  |A|  and  phase  <fr: 

!A!X  =  (-kQi+klr)  |A|  +  Xr|A|3 
*X  =  kli  +  X1  1A|2 

Here  the  subscript  "r"  and  "1”  refers  to  the  real  and  imaginary 
components  of  the  complex  quantities  k  and  X. 

We  shall  briefly  discuss  the  types  of  mathematical  problems  with  which 
we  are  dealing.  In  the  fundamental  mode  equation,  there  are  three 
parameters  appearing  In  the  equations:  R,  w,  kQ.  Since  we  are  Inter¬ 
ested  in  spatially  amplified  solutions  of  the  Orr-Sommerfeld  equation, 
two  of  the  parameters,  R  and  w,  are  given  and  k0  is  to  be  determined. 
In  this  form,  the  fundamental  mode  equation  is  a  nonlinear  problem,  with 
kQ  being  the  eigenvalue  which  is  treated  as  an  additional  dependent 
variable.  The  type  of  system  is  a  two-point  boundary  value  problem 
(2PBVP),  since  the  boundary  conditions  on  (and  its  derivatives)  are 
Imposed  on  both  ends  of  the  computational  Interval. 
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The  adjoint  problem  (2)  has  the  same  form  as  the  fundamental  mode  pro¬ 
blem,  except  that  the  eigenvalue  k0  Is  known  in  advance  from  the  solu- 
tion  of  the  fundamental  mode  problem,  since  the  operators  Lu  and  Lu  can 
be  shown  to  have  the  same  eigenvalue.  On  the  other  hand,  the  second 
harmonic  problem  (5)  is  an  inhomogeneous  2PBVP,  not  an  eigenvalue  prob¬ 
lem.  Finally,  the  correction  to  the  mean  flow  (3),  in  contrast  to  the 
other  three  above-mentioned  problems,  is  a  linear  partial  differential 
equation  (PDE),  similar  in  form  to  the  boundary  layer  equation.  Thus, 
in  addition  to  the  boundary  conditions,  initial  conditions  have  to  be 
specified  for 


The  remaining  two  problems  of  determining  X  and  k^  are  pure  quadra¬ 
tures.  Since  the  integrands  are  going  to  be  approximated  on  some  compu¬ 
tational  grid,  the  integrals  can  be  approximated  by  applying  some  simple 
quadrature  rules,  like  the  trapezoidal  rule,  to  the  integrands  and 
integrating  the  exponential  tails  analytically  from  the  end  of  the  com¬ 
putational  interval  to  infinity. 

The  functions  and  can  be  computed  as  follows.  Note  first 

that  the  differential  operator  L^kg)  appearing  in  problems  (31  and  (7) 

is  singular  and  the  solvability  conditions  simply  force  the  right-hand 

sides  to  be  orthogonal  to  the  eigensolution  of  L  (k„l  so  that  nontrivial 

(11  (11  0 

solutions  exist.  However,  <j>:  and  '  are  only  determined  up  to  the 

1  L  (11  (ll  . 

addition  of  an  arbitrary  multiple  of  $  ,  since  <Jr  '  is  the  eigen¬ 

solution  of  Lu(kol  .  Hence,  we  need  to  specify  normalization  conditions 
for  them  in  order  to  remove  this  degree  of  freedom.  To  be  consistent 
with  the  assumed  form  of  the  expansion  for  the  total  stream  function  ip, 
such  a  condition  could  be  that  and  <$>2^  should  be  orthogonal  to  the 
eigensolution  of  kQ )  : 


(!)*«,.(!) 


B<J>!  dy 


r 0 


0 
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where  8  is  the  operator  defined  by  (see  Appendix  for  the  orthogonality 
condition) 


B  =  4kQ(D2-k2)  +  i 


iR^UfD2- 


3k2)  +  2uk 


,-U"] 


Therefore,  the  desired  solution  can  be  recovered  by 
=  <t>[^( computed)  -  ^ 

(computed)  _  c*X)  *(1) 


where 


c[1}  =  f  (computed)  dy/f  ^l^*B^l'dy 

•^o  '  ■'o 

and  C*”  is  obtained  by  placing  computed)  instead  of  <t>[^  (com¬ 
puted)  . 

From  the  above  brief  description,  we  see  that  at  least  three  of  the 
first  four  problems  are  2PBYP.  It  turns  out  that  if  we  attempt  to  solve 
the  distortion  of  the  mean  flow  equation  by  some  marching  technique 
(e.g.,  the  methods  of  lines,  to  be  discussed  later),  a  sequence  of  2PBVP 
has  to  be  solved.  Since  these  four  problems  constitute  the  major  por¬ 
tion  of  the  overall  computational  problem,  it  becomes  immediately  obvi¬ 
ous  that  we  need  an  accurate  and  efficient  2PBYP  solver  that  is  able  to 
handle  the  different  types  of  2PBVP  that  we  have  discussed  above.  In 
the  next  section,  we  shall  discuss  in  some  detail  the  choice  of  such  a 
solver. 


f 
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3.  NUMERICAL  METHODS  FOR  TWO-POINIT  BOUNDARY  VALUE  PROBLEM 

A  classical  technique  for  solving  the  Orr-Sonwnerfeld  equation  is  to  use 
a  shooting  scheme  together  with  the  Runge-Kutta  method  on  a  uniform  mesh 
(e.g.,  Reynolds  and  Potter  (1967),  Ling  and  Reynolds  (1973)).  The  pre¬ 
sence  of  linearly  independent  solutions  with  vastly  different  growth 
rates,  is  usually  handled  by  some  sort  of  filtering  to  maintain  the  in¬ 
dependence  of  the  computed  solution  (e.g.,  Kaplan  (1964),  Scott  and 
Watts  (1975)).  There  are  also  other  methods  that  do  not  take  the  ini¬ 
tial-value  problem  approach  but  rather  solve  the  boundary  value  problem 
directly  using  finite  difference  methods  (e.g.,  Thomas  (1953),  Osborne 
(1967),  Orszag  (1971)).  Gersting  and  Jankowski  (1972)  discussed  a  com¬ 
parison  of  the  performance  of  these  methods  on  the  Orr-Sommerfeld  pro¬ 
blem. 

The  development  of  modern  numerical  methods  for  general  nonlinear  2PBVP 
has  greatly  advanced  In  the  past  decade  (see,  for  example,  the  surveys 
by  Keller  (1975)  and  Keller  (1976)).  We  now  have  available  powerful 
packages  of  computer  codes  designed  for  a  large  class  of  2PBVP,  of  which 
the  Orr-Sommerfeld  and  related  problems  are  particular  cases.  Many  of 
these  packages  incorporate  adaptive  mesh  selection  techniques  and  can 
handle  rather  sharp  boundary  layers  efficiently.  A  very  high  order  of 
accuracy  can  be  achieved  through  the  use  of  extrapolation  techniques. 
Reliable  error  estimates  are  sometimes  available  and  nonlinearity  is 
routinely  handled.  We  choose  to  use  one  of  these  packages  for  the  2PBVP 
with  which  we  are  concerned. 

The  solver  we  use  is  called  PASVA3,  the  original  version  of  which  is 
described  in  Lentini  and  Pereyra  (1977).  We  shall  briefly  describe  its 
underlying  methods  and  capabilities  here. 
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PASVA3  1$  based  on  the  finite  difference  method  and  Is  applicable  to 
first  order  systems  of  the  form 

=  f(t,y)  ,  t  e  (a,b)  , 

subject  to  the  boundary  conditions  of  the  form: 
gUUhy(b))  =  0  . 

where  f  and  g  can  be  quite  general  nonlinear  functions,  but  sufficiently 
smooth  so  that  the  above  problem  has  an  isolated  solution.  The  method 
of  solution  utilizes  a  global  finite  difference  approximation  rather 
than  a  shooting  method.  Through  the  use  of  the  trapezoidal  rule  approx¬ 
imation,  the  differential  equation  is  approximated  by  a  system  of 
(generally  nonlinear)  algebraic  equations  to  be  solved  for  the  unknown 
function  values  on  the  grid.  Newton's  method  and  variants  of  it  are 
used  to  find  the  solution  to  this  algebraic  system.  The  Jacobian  matrix 
that  results  has  a  special  block  structure  and  an  efficient  elimination 
scheme  is  devised  to  exploit  this  structure.  The  Jacobian  matrix  is 
also  used  to  obtain  a  higher  order  accuracy  (higher  than  the  second 
order  accuracy  of  the  basic  trapezoidal  rule  approximation)  through  the 
use  of  a  procedure  called  deferred  correction.  Since  the  factorization 
of  the  Jacobian  matrix  has  to  be  computed  only  infrequently,  high  order 
accuracies  can  be  achieved  with  very  minimal  costs.  Even  more  important 
than  that  is  the  availability  of  accurate  error  estimates  for  the  com¬ 
puted  solution  through  the  use  of  the  same  Jacobian  matrix.  Accurate 
error  estimates  are  extremely  valuable  because  they  allow  us  to  assign  a 
confidence  level  to  the  computed  results.  Such  error  estimates  are 
noticeably  absent  in  previously  published  computations  related  to  the 
Orr-Sommerfeld  equation.  The  last  significant  feature  of  the  code  that 
we  shall  mention  is  the  ability  to  automatically  and  adaptively  select  a 
good  mesh  for  the  problem.  This  feature  allows  the  code  to  resolve  mild 
boundary  layers  (region  of  relatively  sharp  changes  in  the  solution) 
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efficiently  and  at  the  same  time  reduce  the  computer  storage  required 
for  a  given  accuracy.  A  crude  adaptive  mesh  selection  procedure  was 
first  used  by  Hughes  (1972)  on  the  Orr-Sommerfeld  problem. 


Computational  Domain 


We  used  a  general  rectangular  finite  difference  grid  (Figure  3-1)  on 
which  the  problems  listed  in  Section  1  are  to  be  solved.  This  grid  is 
non-uniform  in  general. 


y, 

'n 

■ 

1 

i 

X 

!  xL 

1  L 

Figure  3-1.  The  Computational  Grid 


On  this  grid,  the  dimensionless  mean  flow  stream  function  T  (and  its  de¬ 
rivatives)  are  to  be  given,  depending  on  the  value  of  the  Reynolds  num¬ 
ber  and  the  kind  of  problem  itself.  Also,  appropriate  initial  condi¬ 
tions  are  to  be  given  at  x  3  x^  and  appropriate  boundary  conditions  to 
be  given  on  y  =  y^  and  y  =  yN.  This  formulation  of  the  computational 
domain  is  quite  general,  and  can  handle  boundary  layer  flows,  Polseuille 
flows  and  Couette  flows.  For  example,  for  plane  Poiseuille  flow  with 
walls  at  yj  =  -1,  yN  =  1,  Tty)  3  y  -  y-  is  known  exactly.  For  Blasius 
flow  with  Reynolds  number  R,  the  dimensionless  mean-flow  stream  func¬ 
tion  T  at  a  point  (x,y)  is  given  by 
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and  f  satisfies  the  following  2PBVP: 

2  f' "  +  ff"  =  0  ) 

>  (3.2) 

f (0)  *  f'(0)  =  0  ,  f'(-)  =  1.  J 

The  function  f(n)  is  calculated  using  PASVA3.  The  mean-flow  on  the 
grid  can  then  be  calculated  from  (3.1)  by  interpolation  from  a  table  of 
f  vs.  n.  The  value  of  yN  is  chosen  so  that  the  corresponding  values  of 
n  at  all  x- stations  are  large  enough  to  contain  most  of  the  boundary 
layer.  Usually,  a  value  of  n  about  10  is  sufficient. 

The  placement  of  the  y  mesh  points  is  determined  by  the  adaptive  mesh 
selection  procedure  in  PASVA3  when  we  use  it  to  solve  the  fundamental¬ 
mode  problem  at  x  =  x1  to  a  given  accuracy.  If  the  mean  flow  does  not 
vary  with  x  too  strongly,  this  mesh  should  be  good  enough  for  all  sub¬ 
sequent  x-statlon  fundamental  modes.  This  mesh  should  also  be  fairly 
good  for  the  second  harmonic  and  the  distortion  of  mean  flow  equations 
because  the  fundamental  mode  constitutes  the  forcing  functions  for  those 
equations.  Our  experience  seems  to  indicate  that  a  mesh  selected  this 
way  is  more  than  adequate  for  all  the  three  other  problems  (Including 
the  adjoint  problem).  Notice  that.  In  principle,  it  is  not  necessary  to 
use  the  same  mesh  in  y  for  all  x-stations;  its  choice  is  for  conveni¬ 
ence,  especially  in  connection  with  the  marching  scheme  for  solving  the 
distortion  of  mean  flow  problem  (6). 
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The  placement  of  the  x-mesh  Is  to  be  specified  by  the  user.  It  can  be 
non-uniform,  but  for  solving  the  problem  (6),  the  x-mesh  should  be  fine 
enough  to  resolve  the  variations  In  the  x-directlon,  of  4>^,k0,k^  and  A 
but  not  too  fine  In  order  to  save  computational  effort. 

The  marching  scheme  for  the  distortion  of  the  mean  flow  equation  will  be 
described  In  Section  3.2  and  the  correct  specification  of  the  boundary 
conditions  will  be  discussed  in  Section  3.3.  In  the  next  section,  we 
shall  describe  how  we  handle  the  eigenvalue  problem  (1). 


3.1  Treatment  of  the  Eigenvalue  Problem 


As  we  have  mentioned,  problem  (1)  is  an  eigenvalue  problem;  both 
and  kQ  are  to  be  determined.  The  approach  taken  by  many  people  is 
to  drop  one  boundary  condition  at  one  end  of  the  computational  interval 
and  impose  a  normalization  boundary  condition  at  the  other  end,  thus 
making  the  new  problem  inhomogeneous,  allowing  a  solution  to  be  easily 
computed  for  a  given  estimate  of  kQ.  An  Iterative  scheme  is  then  used 
to  find  the  correct  value  of  k0  so  that  the  dropped  boundary  condition 
is  satisfied.  The  adjusting  of  k0  can  be  done  by  Newton's  method,  for 
example.  We  have  chosen  another  approach  which,  as  discussed  in  Keller 
(1976),  is  to  view  k0  as  an  extra  unknown  (in  addition  to  <j>)  and  intro¬ 
duce  an  extra  equation 

d 

■  (Jy  '  =  0  (3.3) 

since  k0  Is  constant  Independent  of  y  and  an  extra  normalization  bound- 
(1) 

ary  condition  for  '  such  as 


<yt)  =  i  • 


(3.4) 
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The  "Inflated"  system  becomes  Inhomogeneous  because  of  (3.4)  and  Is  of 
the  form  suitable  for  PASVA3.  In  effect,  both  and  k0  are  obtained 

through  the  use  of  Newton's  method  that  Is  built-in  In  PASVA3.  This 
greatly  simplifies  the  computer  programming  In  the  actual  Implementa¬ 
tion. 


REDUCTION  TO  FIRST  ORDER  SYSTEMS 


PASVA3  accepts  only  real  first  order  systems  of  2PBVP.  However,  the 
problems  listed  in  Section  1  are  all  of  higher  orders  and  are  complex  In 
general.  Therefore,  these  equations  have  to  be  reduced  to  equivalent 
real  first  order  systems  before  PASVA3  can  be  applied.  Except  for  the 
distortion  of  mean  flow  problem  (6),  this  can  be  done  rather  straight¬ 
forwardly.  We  shall  present  these  reductions  In  this  section. 


Fundamental  Mode  Equation 

If  we  define  u  =  4^,  4>y^  .  4>yyy,  kQ)T  ,  then  the  Orr  Sommerfeld 

equation,  together  with  (3.3),  reduce  to  the  following  system  of  five 
first-order  equations: 


du 


u2 

u3 


auj  +  BU3 
0 


=  f(u) 


(3.5) 
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where 


.. .etc. 


further, 


a  =  -  TR(u5  *y  -  «)  u|  -  iR  U5  *yyy  -  ujt 
and 

S  ■  1  R(Ug  1*  -  «)  +  2  u|  . 

To  reduce  (3.5)  to  a  real  system,  we  Introduce  the  real  and  Imaginary 
parts  of  u  as 

JA  *  )L  +  1  l  .  (3.6) 

where  v  and  z  are  real . 

Then,  (3.5)  Is  equivalent  to 

d  m  _[Re  £  (X  +  1  X)] 

w  UJ  l1"1  £  + 1 

where  Re  and  Im  denote  the  real  and 
first  order  system  of  order  10  to  which  PASVA3  can  be  directly  applied. 

The  system  (3.7)  is  nonlinear  and  PASVA3  requires  the  computation  of  the 
Jacobian  of  the  right  hand  side  F  in  (3.7).  From  (3.7)  it  follows  that, 
in  block  form, 


=  F(v,£)  (3-7) 

imaginary  parts.  This  is  a  real 
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J 


3F 

3(v,zT  ' 


a  Re  f  3Re  f 

_32~ 


aim  f  3lm  f 
“Jv-"  “3£— 


t 


► 


» 


I 


Since  fly)  Is  an  analytic  function  In  y  , 
easily  from 

3  f  3£  3  f  3 1 

37"  3i  •  31*  1  3u 


3f 

Now,  ^  Is  easily  computed  from  (3.5)  as 


«y 

'0 

1 

0 

0 

0  ' 

n 

0 

0 

1 

0 

0 

i 

0 

0 

0 

1 

0 

a 

0 

0 

0 

Y 

.0 

0 

0 

0 

0  . 

where  , 

p 

RL*y 

“I’ 

2 

u5lu5*y- 

•0 

R 

+  4 

U5 

]U3 

Hence,  the  Jacobian  of  £  Is  given  by 


each  submatrix  can  be  computed 
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where  the  subscripts  r  and  1  on  a,  6,  and  y  denote  the  real  and  Ima¬ 
ginary  parts,  respectively.  In  the  actual  Implementation,  the  vector 
F  and  the  matrix  J  can  be  easily  calculated  as  a  function  of  y  and  i 
(l.e.,  «),  using  complex  arithmetic  If  that  Is  available. 


Second  Harmonic  Equation 


This  equation  Is  linear  and  Its  reduction  to  a  first  order  system  can  be 
carried  out  In  an  analogous,  but  much  simpler,  way  than  for  the  funda¬ 
mental  mode  equation.  With  *^and  kQ  computed  from  the  fundamental 
mode  problem,  the  forcing  function  F^  on  the  right  hand  side  of  equation 
can  be  easily  calculated.  Now  If  we  define  the  vector 


m  (2)  (2)\t 

y  yy  yyy/ 

where  the  notation  (  denotes  the  tranpose. 


»■  (♦m  *m)  • 


(?\ 

Then  It  Is  easily  seen  that  the  fourth-order  equation  for  '  Is  equlv 
alent  to 
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U2 

u3 

u4 

au^  +  Bu^  + 


where 


and 


a  =  -81  R  Uoya,)  kz  -  21  R  kQ*yyy  -  16  kj 
e  =  21  R  (k0* -»)  +  8  kj 


(3.9) 


(3.10) 

(3.11) 


To  further  reduce  (3.9)  to  a  real  system,  we  again  define  the  real  and 
Imaginary  parts  of  u  as  y  *  v  +  1  J  and  rewrite  (3.9)  as 


d 

W 


v2 

v3 

v4 

Re( au^+BU2p2 ) 
z2 
z3 
z4 

Im( Bu^ +8U2+F2) 


The  .Jacobian  matrix  is  easily  computed  from  (3.9)  as 
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J  = 


01000000 
00100000 
0  0  0  1  0  0  0  0 

ar  0  8r  0  -a.  0  -0f  0 


0  0  0  0 
0  0  0  0 
0  0  0  0 
0  6i  0 


0  10  0 

0  0  10 

0  0  0  0 

<v  0  er  0 


(3.12) 


where  the  subscripts  r  and  1  denote  the  real  and  imaginary  parts,  re¬ 
spectively,  of  a  and  3  as  defined  in  equation  (3.10)  and  (3.11). 


Adjoint  Equation 


Analogous  to  the  previous  two  equations,  the  reduction  of  the  adjoint 
equation  can  be  carried  out  as  follows.  Defining 


u  = 


we  get 


W 


u2 

U3 

u4 

aUj+0U2+YU.j 


where  a  =  -i  R(kQijjy-u>)  kjjj-kg 

a  =  2  i  kQ^yy 
y  =  i  R  (kQ^y-w)  +  2  kj 


) 

i 


and 


(3.13) 
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Defining  )£  =  v  +  i  £  as  before,  the  real  system  is 

~  v2 
v3 

Fv1  v4 

W  z  =  Re(au1+Bu2nu3) 

z3 

z4 

Im{au1+Su2+Yu3) 

and  the  Jacobian  matrix  is 


0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

a 

e„ 

Y„ 

0 

-a. 

-8i 

1 

O 

1 

0 

J  = 

r 

0 

r 

0 

r 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

a . 

1 

8i 

Yi 

0 

a 

Y 

0 

r 

r 

r 

using  the  same  notation  as  before. 


(3.14) 


(3.15) 


Notice  that  kQ  Is  a  known  quantity  (obtained  from  the  fundamental  mode 
computation)  In  the  above  formulation,  and  thus  the  problem  is  linear. 
To  make  the  problem  Inhomogeneous,  a  boundary  condition  at  one  end  of 


stm  *  ******* 
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the  computational  interval  is  dropped,  and  an  extra  inhomogeneous  nor¬ 
malization  boundary  condition  is  imposed  at  the  other  end  of  the  inter¬ 
val.  Since  the  fundamental  and  the  adjoint  problems  have  the  same 
eigenvalue,  the  dropped  boundary  condition  should  be  automatically 
satisfied.  This  condition  should  be  checked  after  the  computation. 
More  details  on  boundary  conditions  will  be  discussed  in  Section  3.3. 


3.2  Distortion  of  Mean  Flow  Equation 


This  problem  is  an  initial  and  boundary  value  problem  (IBVP)  for  a  par¬ 
tial  (rather  than  ordinary)  differential  equation.  Since  the  equation 
is  of  parabolic  type  and  the  variation  in  the  x-direction  is  not  expec¬ 
ted  to  be  large,  implicit  schemes  (e.g.,  Crank-Nicol son)  are  more  effi¬ 
cient  than  explicit  methods.  Implicit  methods,  however,  necessitate  the 
solution  of  a  linear  algebraic  system  of  equations  at  each  x-station. 
PASVA3  can  be  used  to  solve  this  linear  system,  thus  utilizing  the 
service  of  its  efficient  elimination  scheme.  Even  more  important  than 
this  is  the  achievement  of  high  order  accuracies  through  the  deferred 
correction  procedure  built-in  in  PASVA3.  We  shall  next  show  how  this 
method  can  be  applied  to  the  distortion  of  mean-flow  equation  and  how 
PASVA3  can  be  used  in  conjunction  with  the  method. 


First  notice  that  the  right  hand  side  Fj  in  the  problem  (6)  is  real  by 
definition  and  so  is  .  We  first  approximate  the  x-derivatives 


of  <(>^  by 


two-point  centered-difference  approximations,  centered 
also  approximate  the  coefficients 


at  xn+l/2=  (xn  *  xn.l)/2  and 

at  xn,.  „  by  averaging.  Specifically,  we  approximate  the  partial  dif- 
n+i/^  / q \ 

ferential  equation  for  <{>  by 
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,n+l 


[C1  »(0)"T‘  *  kSt  ♦(°’n) 


knrl  (0)ntl  +  k 

Koi  vy  9y  K 


n.  /  4>{0)nN) 
oi  ry  y  j 


li)  (,yOI"tl-  ^y0'")"  K^‘ 


1 

JR 


where  the  superscripts  n  and  n+1  correspond  to  the  x-stations.  Now  mul¬ 
tiplying  by  2R  and  collecting  terms  with  superscript  (n+1)  on  one  side, 
we  get 


i>(0) 

yyy 


n+1 


+a 


n+1  (0) 
yy 


(3.16) 


where 


H 


n+1 


♦«»"-  a?  *{0) 
yyy  3  Yyy 


n 


-  R  [Fj+1  + 


F1 1 


(3.17) 


where  the  coefficients  in  (3.16)  and  (3.17)  are  given  by 
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(3.18) 


These  coefficients  are  all  known  quantities  since  If  and  its  derivatives 
are  given  quantities.  The  equation  (3.16)  (with  appropriate  boundary 
conditions)  is  in  the  form  of  a  linear  2PBVP,  with  the  right  hand  side 
function  Hn+1  known  at  the  y-grid  points  at  x  =  xn.  Thus,  PASVA3  can  be 
called  on  to  solve  this  2PBVP  if  we  keep  the  same  y-grid  at  x  =  xn+1  as 
at  x  =  xn.  Of  course,  In  this  fashion,  we  cannot  exploit  the  adaptive 
mesh  selection  capability  of  PASVA3,  unless  we  provide  some  accurate 
interpolation  routine  for  Hn+1.  As  explained  before,  using  the  same  y- 
grid  for  all  x-statlons  Is  adequate  for  near-parallel  flows  and  the  com¬ 
puter  programmi ng  is  kept  simple  this  way.  On  the  other  hand,  high 
order  accuracies  (in  y)  can  still  be  achieved  by  the  deferred  correction 
procedure  in  PASVA3,  though  the  accuracy  in  the  x-discretization  is 
second  order. 

A  reduction  of  (3.16)  to  first  order  system  is  necessary  before  PASVA3 
can  be  applied.  If  we  define 

,  .  r.(o)n+i  .(ointi 

l  -  ?  »  ?>y  » 


i 

i 


then  (3.16)  reduces  to 
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dz 


L2 

z3 


_  -c1  -(*r*  -r1) 


and  the  Jacobian  matrix  is: 


,2  -  af  Zj  a  H"*1 


(3.19) 


J  = 


0 

0 


1 

0 


.n+1 


-(•TX*1) 


a 


0 

1 

n+1 


(3.20) 


3.3  Boundary  Conditions 


The  specification  of  the  problems  will  not  be  complete  without  appro¬ 
priate  boundary  and  initial  conditions.  These  conditions  depend  on  the 
type  of  mean  flow  and  the  computational  domain  of  a  specific  problem. 
In  this  section,  we  shall  consider  two  kinds  of  flows:  plane  Poiseuille 
flow  and  Blasius  flow  over  a  flat  plate. 


3.3.1  Poiseuille  Flow 


Let  the  walls  be  at  y  5  ±1.  Then  the  computational  domain  could  be 
either  [-1,1]  or  the  half-width  [0,1]  for  solutions  symmetric  with  re¬ 
spect  to  the  mid-channel.  The  no- si  ip  condition  at  the  walls  reduces  to 
the  condition  that  the  stream  function  and  its  first  derivative  should 
be  zero  there. 

For  the  domain  [0,1],  the  appropriate  boundary  conditions  are,  for  the 
sywnetrlc  fundamental  mode  problem: 


-37- 


*{y n(°)  =  4yy(0)  3  0  ) 

/  (3.21) 

*(1)(i)  =  *Jn(i)  =  0  ) 

and  the  normalization  boundary  condition  could  be  taken  to  be 

*{1)(0)  =  1  (3.22) 


Then  the  second  harmonic  mode  is  anti -symmetric  about  y  =  0,  and  the 
appropriate  boundary  conditions  are 


*(2)(0)  =  *iy)(0)  =  0 
*(2)(1)  ■  ^1}(1)  *  0  • 


(3.23) 

(3.24) 


The  boundary  conditions  for  the  adjoint  problem  are  the  same  as  for  the 
fundamental  mode  problem. 


For  the  distortion  of  mean  flow,  the  appropriate  boundary  conditions 
on  should  be: 


>(0) (0)  =  ^(0)  *  0 


k(0) 


(3.25) 


(1)  =  0  . 


For  the  perturbed  mean  flow  in  a  channel,  we  can  specify  either  (i)  no 
change  in  the  total  volume  flow  through  the  channel  or  (ii)  no  change 
in  the  pressure  drop  between  two  stations,  x  =  x^  and  x2  say,  along  the 
channel.  For  the  case  (1)  the  boundary  condition  is 


4>(0)  (x,l)  =  0 


(3.26) 
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and  then  the  partial  differential  equation  for  contains  a  non-zero 
function  P(x)  (negative  pressure  gradient)  whose  value  Is  to  be  deter¬ 
mined  so  that  the  last  mentioned  condition  for  Is  satisfied.  For 
the  case  (ii)  the  boundary  condition  for  <t>^  Is 

4>(0){x,l)  =  ^  ,  Independent  of  x  (3.27) 


and  the  value  of  P(x)  In  the  PDE  Is  determined  to  satisfy  the  above 
boundary  condition  everywhere  along  the  channel.  The  value  of  ^  Is  to 
be  determined  so  that  the  imposed  pressure  condition  Is  satisfied: 


f 

X1 


P(x)dx  =  0 


(3.28) 


A 

Only  In  the  temporal  problem,  in  which  the  mean-flow  perturbation  is 
Independent  of  x,  the  last  condition  reduces  to 


P(x)  =  0 


(3.29) 


3.3.2  Blasius  Flow 


At  the  wall,  the  no-slip  boundary  conditions  become 

4>(0)  -  *'(0)  =  0  (3.30) 

where  +  Is  any  of  the  functions  4>(2),  and  .  However, 

the  boundary  conditions  at  the  other  end  of  the  computational  interval 
(l.e.,  y  =  yN),  which  Is  physically  of  the  form 

11m  <#>(y )  =  0  ,  (3.31) 


requires  more  careful  treatment. 
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First  of  all,  the  computational  interval  should  really  be  semi- infinite, 
but  for  practical  reasons,  a  finite  interval  must  be  used  instead. 
Second,  the  perturbations  vanish  at  Infinity.  The  fundamental  mode 
equation,  when  evaluated  outside  the  boundary  layer  has  four  linearly 
Independent  basic  solutions,  two  of  which  decay  while  the  other  two 
grow,  all  exponentially.  We  want  our  computed  solution  at  the  finite 
end  point  to  be  composed  of  only  the  two  decaying  basic  solutions  in 
order  to  satisfy  condition  (3.31).  This  "outgoing  wave"  condition  can 
be  derived  quite  easily  for  the  Orr-Sommerfeld  problem.  The  following 
brief  description  is  adapted  from  Keller  (1976). 

Consider  a  nth-order  linear  2PBVP  of  the  form 
du 

^=A(y)u  (3.32) 

on  the  semi-infinite  Interval  y  >_  0  ,  where  the  limit 

lim  A(y)  =  A*,  (3.33) 

y-H» 

is  assumed  to  exist. 

Let  the  eigenvalues  of  A.  be  denoted  by  X^,  let  and  be  the  inde¬ 
pendent  left  and  right  eigenvectors  of  A*  corresponding  to  x^  ;  that  Is, 

(3.34) 


Then 

(X.-X.)  *T  rt  =  *1  A  r,  -  *T  A  r,  =  0 

1  j  ~j  H  ~j  00  ~1  ~j  »  M 


£I  =  X.  l]  , 


A  r •  =  A  *  fj  « 
00  j  -i 


and  hence 
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t 


# 


t 


> 


r.  =  0  if  1  *  j  , 
~j  ~i 


*  0  if  i  =  j  . 


(3.35) 


Now  every  solution  of  djj/dy  =  Aji  has  the  form 
n 


«(y)  =  ]C  a<Ci 
j=l  3  3 


.V 


(3.36) 


To  ensure  that  the  solution  decays  at  °°,  we  have  to  require  aj  =  0  for 
those  j's  with  Re  A.  0  .  By  the  use  of  (3.35),  this  condition  can  be 
expressed  as: 


Tiro  y(y)  =  0  ,  If  Re  A.  _>  0  . 

y+°o  J  ^ 

Therefore,  the  proper  boundary  condition  for  the  finite  problem,  with 
right  end  point  b,  should  be 


l]  u(b)  =0  if  Re  A.  >  0  (3.37) 

~j  ~  j  — 

Note  that  the  often-used  boundary  condition  of  setting  some  components 
of  u(b)  to  zero  is  inaccurate  if  the  desired  solution  decays  very 
slowly,  unless  the  value  for  b  Is  taken  to  be  large  enough,  which  is 
obviously  undesirable  for  efficiency  considerations. 


For  the  fundamental  mode  problem,  if  we  define 


y  = 


t 


then  we  have 


dy 

<3y  =  AW  M  • 


where 


» 
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A(y)  = 


0 

0 

0 

a 


1 

0 

0 

0 


0  0 
1  0 
0  1 
0  0 


and 


Since 


o  =  -  i  R  (k  tfi  -  u)  -  i  R  k  <ii  -  k* 
*  '  o  vy  '  o  "  *o  *yyy  o 


3  ■  *  R  <kQ  ?  -  «)  +  2  k£  . 


lim  ♦  (y)  *  1  and  Tim  ♦  (y)  =  0 

y+oo  +  y-f«o  JJ'J 


for  Blasius  flow. 


A_  = 


0  10  0 

0  0  10 

0  0  0  1 

a  0  0_  0 


where 


a  =  -  i  R  (k  -w)  k2  -  k^ 
•  0  0  0 


*  <  V1  +  2  ko 


The  eigenvalues  of  K,  are  easily  calculated  to  be 


=  ^9  =  A,  =  Y»  ^4  =  ~T 


(3.38) 


y  =  y  k2  +  i  R  (ICq-m)  »  and  Re(k0)»  Re(T)  >  0 


where 
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The  corresponding  left  eigenvectors  are 


M 

k  / 
V 

—  — 

-k2Y 

V 

V 

-Y2 

-Y2 

-k2 

-k2 

~1 = 

ko 

»~2 

"ko 

*~3 

Ko 

Y 

*~4  = 

Ko 

-Y 

1 

1 

1 

1 

The  eigenvalues  ^  and  have  a  positive  real  part.  The  corresponding 
left  eigenvectors  should  be  used  in  the  boundary  condition  (3.37). 


In  terms  of  the  original  dependent  variable,  and  with  D  =  d/dy, 
condition  (3.37)  reduces  to 


(D2-  k2)  (D  +  y)  =  0 
(D2  -  Y2)  (D  +  kQ)  4>{1)  =  0 


at  y  =  yN  . 


(3.40) 


As  Keller  pointed  out,  these  boundary  conditions  seem  to  have  rarely 
been  used  other  than  in  shooting  methods  for  the  Orr-Sommerfeld  prob¬ 
lems,  although  related  form  such  as 


(D  +  y)  =  0 
(D  +  kQ)  =  0 


at  y  =  yN 


at  y  =  yN 


(3.41) 


(3.42) 


have  been  used.  Grosch  and  Orszag  (1977)  have  found  that  conditions 
like  (3.41)  and  (3.42)  can  perform  rather  poorly. 
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In  our  formulation  where  kQ  Is  treated  as  an  extra  unknown  as  shown  In 
Section  5.1,  an  extra  normalization  boundary  condition  like 

♦^(0)  =  1 

have  to  be  Imposed.  This  completes  the  description  of  the  boundary  con¬ 
ditions  for  the  fundamental  mode  problem. 

For  the  adjoint  problem.  It  can  be  easily  verified  that  the  matrix  A.  Is 
the  same  as  for  the  fundamental  mode  problem  (l.e.,  the  Orr-Sommerfel d 
problem  Is  self-adjoint  at  ®) ,  and  hence  the  same  boundary  conditions 
(3.40)  should  be  used. 

Boundary  conditions  for  the  second  harmonic  equation  and  the  distortion 
of  mean  flow  equation  can  be  derived  In  a  similar  fashion.  The  presence 
of  forcing  terms  In  the  differential  equation  requires  a  slightly  dif¬ 
ferent  treatment. 
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4 .  NUMERICAL  RESULTS 

The  numerical  procedure  described  in  the  previous  sections  can  be  ap¬ 
plied  to  very  general  classes  of  hydrodynamic  stability  problems.  Our 
main  Interest  In  the  present  report  is  to  compute  the  solution  for  a 
particular  case  of  boundary-layer  stability  over  a  flat  plate,  and  com¬ 
pare  the  numerical  results  to  those  computed  by  Murdock  (1977)  using  a 
totally  different  approach.  However,  in  order  to  have  some  feeling  for 
the  performance  and  accuracy  of  the  present  numerical  procedure,  we  also 
applied  the  procedure  to  two  cases  where  published  results  are  avail¬ 
able. 

4.1  Plane  Polseullle  Flow  (Reynolds  and  Potter  (1967)) 

Computations  were  performed  for  the  case  of  R  =  5772.12,  u>  *  0.2694879.* 
This  corresponds  to  the  critical  point  on  the  neutral  curve.  Reynolds 
and  Potter  (1967)  also  computed  solutions  for  this  case,  but  from  a 
temporal  stability  standpoint.  In  which  R  and  kQ  are  given  and  * 
(complex,  in  general)  Is  computed  as  the  eigenvalue  of  the  Orr- 
Sommerfeld  problem.  But  for  points  on  the  neutral  curve,  w  Is  purely 
real  so  the  temporal  solution  Is  identical  to  the  spatial  stability  so¬ 
lution  obtained  by  specifying  a  real  value  for  u>  and  computing  the 
generally  complex  k0  as  the  eigenvalue.  The  wavenumber  kQ  used  by 
Reynolds  and  Potter  Is  1.02071  and  the  frequency  w  as  computed  by  them 
Is  0.2694879  +  01.  Therefore,  we  solved  the  case  with  «  «=  0.2694379 
(real)  as  Input,  and  the  eigenvalue  kQ  as  computed  by  us  Is 

1.0207099  +  4.6032093  x  10'7  1  , 

which  agrees  with  Reynolds  and  Potter's  result  to  all  significant  digits 
given. 


differs  from  theirs  by  a  factor  of  3/2. 
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Pgiseuille,  R=5772.12,6j=.26949 
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Poiseuille,  R=5772.12,u=.26949 


o 
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This  particular  computation  was  done  on  a  half-depth  domain  of  [0,1] 
starting  with  50  uniformly  spaced  grid  points.  PASVA3  added  31  more 
points  and  took  about  2  seconds*  to  achieve  an  estimated  relative  error 
of  10"^  in  4^^  and  all  its  derivatives  up  to  the  third.  As  discussed 
before,  PASVA3  gives  an  estimate  for  the  error 


max  1  u  -  ( computed)  -  u exact)  1  ,  y  Ety^y^] 

for  each  component  of  the  computed  solution  u. 
estimated  relative  error  is  defined  as 


In  this  report,  the 


max  |u.j(computed  -  u^(exact)! 
max  fiiT (computed) I  ‘ 

The  initial  guess  for  the  and  its  derivatives  are  simply  constants 
and  not  close  to  the  exact  solutions  at  all.  The  initial  guess  for  kQ 
was  unity.  PASVA3  encountered  no  difficulties  with  convergence  to  the 
correct  eigenvalues  and  eigenfunctions.  The  computed  ,  and 

<t>^  are  plotted  in  Figures  (4.3)-(4.3)  where  we  show  only  every  other 
grid  point.  It  is  seen  that  PASVA3  does  indeed  put  more  grid  points  in 
places  where  the  functions  change  the  most.  These  plots  are  indistin¬ 
guishable  from  those  computed  by  Reynolds  and  Potter. 


4.2  Linear  Stability  of  Blasius  Flow  (Jordinson  (1970)) 


We  also  computed  the  stability  solutions  for  Blasius  flow  over  a  flat 
plate  and  compared  them  to  three  cases  computed  by  Jordinson. 
Jordinson' s  three  cases,  after  transformation  to  our  dimensionless  vari¬ 
ables  (in  our  computation  we  scale  lengths  by  the  distance  from  the 
leading  edge  whereas  Jordinson  scales  by  the  displacement  thickness), 
are  tabulated  in  Table-4.1. 
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Table  4.1. 


Case 

R 

CO 

kor 

koi 

1 

38125.69 

14.717 

34.994 

0.8964 

2 

120765.14 

24.254 

62.1799 

-0.3837 

3 

336356.86 

37.815 

104.0077 

-1.9211 

Our  computed  results  are  tabulated  in  Table  4.2. 


Table  4.2 


Case 

kor 

koi 

ymax 

nmax 

N 

Approximate 

Relative  Error 

R  R  • 

'•or  Koi 

1 

34.98617 

0.905057 

0.05 

9.76 

99 

1 

O 

rH 

10-4 

2 

62.17120 

-0.374498 

0.03 

10.43 

95 

10'6 

1 

o 

3 

104.01196 

-1.907428 

0.015 

• 

00 

97 

10-6 

1 

o 

All  computations  were  started  with  20  uniformly  spaced  grid  points  and 
ymax  is  chosen  large  enough  so  that  it  is  effectively  outside  the  bound¬ 
ary  layer  thickness.  The  execution  time  for  each  case  was  about  1.5 
seconds.  Jordinson  used  a  finite  difference  technique  with  80  uniformly 
spaced  grid  points  and  he  claimed  an  accuracy  of  about  five  decimal 
places  which  is  equivalent  to  a  relative  error  of  about  10“5  for  the 
three  cases  above.  But  on  comparison  of  the  results  shown  in  Tables  4.1 
and  4.2,  we  see  that  our  results,  in  general,  agree  with  Jordinson' s  to 
about  the  same  accuracy  for  kQr  but  our  k0j's  differ  by  much  more. 
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Direct  comparison  with  Jordinson's  eigenfunctions  is  not  possible  be¬ 
cause  we  do  not  know  how  he  normalizes  his  eigenfunction.  Without  the 
availability  of  an  exact  solution,  it  is  not  possible  to  conclude  which 
computation  is  more  accurate.  However,  the  availability  of  reliable 
error  estimates  in  our  computation  certainly  gives  us  more  confidence  in 
the  estimated  accuracy  of  our  computed  results,  especially  when  viewed 
in  conjunction  with  the  accuracy  of  our  Poiseuille  flow  calculations  in 
Section  4.1.  It  is  precisely  in  such  situations  when  exact  solutions 
are  not  easily  available  that  reliable  error  estimates  like  those  com¬ 
puted  by  PASVA3  are  invaluable  in  assessing  the  accuracy  of  the  computed 
resul ts. 


4.3  Nonlinear  Stability  of  Blasius  Flow  (Murdock  (1977)) 

The  main  goal  of  our  present  effort  is  to  compute  the  growth  of  distur¬ 
bances  as  they  propagate  downstream  from  the  stable  region  through  the 
neutral  point  to  the  unstable  region.  Murdock  (1977)  computed  the 
growth  of  a  large  amplitude  Tollmien-Schl ichting  wave  for  the  case: 
R  =  105,  u  =  13.19  from  x  =  1.0  to  x  =  2.2  in  our  dimensionless  vari¬ 
ables.  His  approach  was  to  solve  the  time-dependent  "parabolized" 
Navier-Stokes  equation  with  the  Tollmien-Schl ichting  wave  as  initial 
condition.  In  order  to  compare  with  his  results,  we  computed  for  this 
case  the  quantities  <|>(1  ,4>(0)  .♦J1J  ,kQX,ki  ,  and  ^  , 

on  a  finite  difference  grid  with  96  non-uniform  mesh  points  (N  =  96), 
13  x-stations  (L  =  13)  uniformly  spaced  between  x  =  1.0  and  3.4.  yN  w.is 
chosen  to  be  0.04  so  that  ^  =  12.6.  The  mesh  in  y  was  determined  by 

solving  the  Orr-Sommerfeld  problem  with  initially  20  uniformly  spaced 
grid  points  in  the  y-direction.  The  estimated  relative  error  in  the 
numerical  solutions  as  computed  by  PASVA3  are  all  less  than  10"3  and,  in 
most  cases,  less  than  10‘4.  To  give  an  idea  of  the  efficiency,  the  exe¬ 
cution  times  are  summarized  in  Table  4.3. 


i 

CPU  seconds  per  x-station  j 

*U> 

nonlinear 

< 

1 

♦<2) 

1  inear 

< 

.25 

*(1) 

X 

linear 

< 

.25 

*i 

linear 

< 

.25 

A 1) 

*2 

1  inear 

< 

.25 

AD* 

9 

1 i near 

«S 

.3 

JO) 

4> 

linear 

< 

.1 

A,k1 

< 

.1 

Table  4.3 


The  results  of  our  computation  for  uj  =  13.19  are  summarized  in 
Table  4.4.  X  and  X.  vary  somewhat  irregularly  at  a  few  points  near 
x  =  1  and  this  behavior  probably  was  caused  by  our  choice  of  <j>  '  =  0 

(correction  to  the  mean  flow)  at  x  =  0*. 


c 

RxlO  D 

kor 

koi 

1.0 

36.67 

1.079 

1.2 

36.32 

.398 

1.4 

36.08 

-.059 

1.6 

35.91 

-.360 

1.8 

35.80 

-.545 

2.0 

35.71 

-.642 

2.2 

35.65 

-.668 

2.4 

35.60 

-.639 

2.6 

35.56 

-.562 

2.8 

35.53 

-.447 

3.0 

35.49 

-.298 

3.2 

35.45 

-.120 

3.4 

35.41 

+  .084 

xrxl03 


x.xlO3 


-.422 

-.349 

-.294 

-.252 

-.218 

-.191 

-.168 

-.148 

-.131 

-.115 

-.101 


-1.22 

-1.26 

-  .97 

-  .67 

-  .49 

-  .20 

+  .15 

.73 

1.56 

2.76 

4.29 

6.12 

8.02 


-  .15 
-1.70 
-1.24 
-1.84 
-2.05 
-2.67 
-3.27 
-4.19 
-5.17 
-6.34 
-7.44 
-8.44 
-9.14 


Table  4.4  Stability  Parameters  (< 


We  made  this  choice  in  order  to  compare  our  computation  with  Murdock's 
(see  next  section) . 
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For  infinitesimal  disturbances,  kor  and  k0j  are  the  wave  number  and  the 
damping  rate  in  local  parallel  mean  flow.  kjr  and  k^-  represent  the 
corrections  to  the  wave  number  and  amplification  rate  resulting  from  the 
slow  variation  of  the  mean  flow  with  the  streamwise  distance.  If  we 
define  the  corrected  amplification  rate  by  -kQi  +  klr  as  in  Saric  and 
Nayfeh  (1975),  the  neutral  Reynolds  number  (corresponding  to 
-k01-  +  klr  =  0)  is  decreased  by  5.4%  from  that  for  the  local  parallel 
flow  (kQ1-  =  0),  in  close  agreement  with  Saric-Nayfeh' s  result. 

x  (=x  +ix.)  Is  the  Landau  constant  which  is  the  measure  of  nonlinear 

r  1  (2) 
effects  on  the  stability.  The  second-harmonic  component,  r  ,  and  the 

correction  to  the  mean  flow,  contribute  nearly  equally  to  the  con¬ 

stant  x. 


With  the  values  of  kQ,  kj,  X  computed  for  a  range  of  x,  the  "amplitude" 
A(x)  of  the  basic  component  may  be  computed  from  the  differential  equa- 
ti  on 

^  =  [i  kQ(x)  +  kx(x)]  A  +  x ( x )  I A 1 2  A 


The  solution  is  given  by 
|A( 


x) I  =  AQ  p(x)  j jl  -  A2  q(x)j 


1/2 


A(x)  =  1 A( x 


)!  exp  |  i  j  (kQr  +  kXi  +  x.|A|2)  d£ 
xo 

■I  ft 


p(x)  =  exp<|  |  |-k1„(C)  +  k,„U)|  dC 


r 


q(x)  =  2  J  Xr(C)  p2U)  dK  . 


|A(xq)|  . 


where 
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Thtis,  If  A.^ I q( x)  | <<1 ,  we  get 

i-  I A{  x )  |  =  p( x ) 

Ao 

which  is  the  linear  result  discussed  earlier.  The  nonlinear  effects 
reduces  the  amplitude  in  the  region  where  q(x)<0,  and  Increases  the  am¬ 
plitude  where  q( x) >0-  When  A*  q{x)  becomes  unity,  |A(x)|  becomes  in¬ 
finite.  A  similar  instability  "burst"  was  reported  by  Hocking,  Stewart- 
son  and  Stuart  (1972)  in  their  study  of  infinitesimal  three-dimensional 
wave  packet  in  a  fully  developed  plane  Poiseuille  flow  at  a  slightly 
supercritical  Reynolds  number.  In  a  slowly  changing  boundary-layer 
flow,  q(x)  is  also  a  slowly  changing  function  of  x  which  is  negative  for 
low  Reynolds  number  (based  on  x)  region  and  becomes  positive  beyond  a 
certain  Reynolds  number.  It,  however,  is  a  universal  function  for  a 
fixed  initial  station  xQ  in  a  given  undisturbed  boundary  layer,  and  once 
computed  it  may  be  used  to  estimate  the  "burst"  location  as  AQ  is 
varied. 

The  amplitude  functions  p(x)  and  q(x)  for  w  =  13.19  are  given  in 
Table  4.5.  The  last  column  gives  the  value  of  A0  for  which  the  corres¬ 
ponding  x  is  the  singularity  of  A(x).  The  numerical  value  for  A  depends 
on  how  the  eigenfunctions  are  normalized.  In  our  computation  we  arbi¬ 
trarily  chose  <J>'  '(0)  =  1,  and  I  <P'  (y)  !max  (corresponding  to  AQ  =  1  at 
x  =  1)  is  about  0.15  x  10”  2 .  Hence,  though  A0  =  0(10)  seems  large,  the 
streamwise  velocity  fluctuation  is  only  0(3%).  Near  the  burst,  the  am¬ 
plitude  becomes  so  large  and  the  weakly  nonlinear  theories  become  in¬ 
valid.  If  we  assume  that  truly  nonlinear  processes  leading  to  boundary 
layer  transition  take  place  in  a  short  distance  around  the  Instability 
burst  predicted  by  the  present  theory,  then  we  can  use  the  burst  loca¬ 
tion  as  an  estimate  for  the  transition  location  replacing  the  en-cri- 
terion.  For  this  purpose,  we  have  to  compute  q(x)  for  a  range  of  fre^ 


*  This  instability  "burst"  should  not  be  confused  with  the  instability 
burst  observed  experimentally  in  boundary  layers.  The  "burst"  in  the 
present  case  simply  means  the  break-down  of  weakly  nonlinear  solu¬ 
tion. 
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X 

p(x) 

q(x)*100 

Ao 

1. 

1. 

0 

★ 

1.05 

.96 

-.012 

★ 

1.1 

.93 

-.023 

* 

1.15 

.91 

-.034 

* 

1.2 

.89 

-.044 

★ 

1.25 

.88 

-.054 

★ 

1.3 

.88 

-.063 

★ 

1.35 

.88 

-.072 

★ 

1.4 

.89 

-.08 

★ 

1.45 

.9 

-.087 

★ 

1.5 

.92 

-.094 

★ 

1.55 

.94 

-.101 

★ 

1.6 

.96 

-.107 

* 

1.65 

.99 

-.113 

★ 

1.7 

1.02 

-.119 

★ 

1.75 

1.05 

-.125 

★ 

1.8 

1.09 

-.131 

★ 

1.85 

1.13 

-.137 

★ 

1.9 

1.17 

-.142 

★ 

1.95 

1.22 

-.147 

★ 

2. 

1.27 

-.151 

★ 

2.05 

1.32 

-.153 

★ 

2.1 

1.38 

-.155 

★ 

2.15 

1.44 

-.154 

★ 

2.2 

1.5 

-.152 

★ 

2.25 

1.56 

-.147 

★ 

2.3 

1.63 

-.138 

★ 

2.35 

1.69 

-.125 

★ 

2.4 

1.76 

-.106 

★ 

2.45 

1.83 

-.079 

★ 

2.5 

1.91 

-.044 

★ 

2.55 

1.98 

2E-03 

212.7 

2.6 

2.05 

.061 

40.6 

2.65 

2.13 

.134 

27.27 

2.7 

2.2 

.227 

21.01 

2.75 

2.27 

.34 

17.15 

2.8 

2.34 

.478 

14.46 

2.85 

2.41 

.644 

12.46 

2.9 

2.48 

.842 

10.9 

2.95 

2.55 

1.075 

9.65 

3. 

2.61 

1.346 

8.62 

3.05 

2.66 

1.659 

7.76 

3.1 

2.72 

2.018 

7.04 

3.15 

2.76 

2.425 

6.42 

3.2 

2.81 

2.882 

5.89 

3.25 

2.84 

3.389 

5.43 

3.3 

2.87 

3.946 

5.03 

3.35 

2.89 

4.553 

4.69 

3.4 

2.91 

5.206 

4.38 

Table  4.5.  Nonlinear  Stability  Amplitude  Functions 


Since  the  consensus  in  the  transition  research  community  Is  that  three- 
dimensional  disturbances  are  more  important  than  two-dimensional  distur¬ 
bances  In  pre-transition  nonlinear  regions,  a  more  pressing  task  is  to 
extend  the  present  theory  to  weakly  nonlinear  three-dimensional  distur¬ 
bances. 


Comparison  with  Murdock's  Computation 

Murdock  (1977)  investigated  the  nonlinear  effects  on  the  stability  of 
boundary  layer  over  a  flat  plate,  by  numerically  integrating  the  non¬ 
steady  Navier-Stokes  equations.  The  computations  are  carried  out  in  the 
Reynolds  number  range  10s  <_  U^x/v  £  2.5  x  105  .  At  the  upstream  bound¬ 
ary!  U^x/v  =  105,  the  condition  Is  the  superposition  of  Blasius  steady 
boundary-layer  solution  and  a  time-periodic  solution  of  the  temporal 
Orr-Sommerfeld  equation.  He  presents  the  results  for  dimensionless  fre¬ 
quency  w  Xj/U^  =  13.19  and  amplitudes*  A  =  0.001  and  0.08. 

Our  result,  corresponding  to  Murdock's  A  -  0.001,  Is  shown  in 
Figures  4.4a,b.  We  determined  the  initial  amplitude  by  matching  his 
results  at  the  initial  station.  Since  he  used  the  temporal  solution  of 
the  Orr-Sommerfeld  equation  at  the  upstream  boundary  while  our  solution 
is  the  spatial  solution  (the  frequency  w  is  real  and  the  wave  number  k 
is  complex),  a  perfect  match  was  not  possible  and  we  used  the  arithmetic 
average  of  two  values  obtained  by  matching  at  y^U^72vx  =  0.2  and  1.  In 
the  case  corresponding  to  Murdock's  A  =  0.001,  the  nonlinear  effects  are 
negligible  and  the  amplitude  is  given  by  the  solution  of 


1 

m 


d  1 A I  _ 


'koi  +  klr 


*  The  meaning  of  A  is  not  clearly  defined  in  the  paper,  since  he  failed 
to  state  how  the  Orr-Sommerfeld  solution  is  normalized. 
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f  Figure  4.4a  shows  the  variation  of  the  amplitude  of  streamwise  component 

of  perturbation  velocity  along  y  =  0.2-^2v  x/U^.  The  agreement  between 
our  result  and  Murdock's  computation  Is  only  fair,  and  we  feel  that  the 
difference  is  caused  by  the  imperfect  matching  of  the  boundary  condition 
at  x  =  0. 

The  comparison  in  the  large  amplitude  case  (A  =  0.08  in  Murdock)  turned 
out  to  be  more  complicated.  Murdock  simply  increased  the  amplitude  of 
the  Orr-Sommerfeld  solution  without  introducing  any  nonlinear  contribu¬ 
tions  at  x  =  1,  whereas,  in  our  formulation  presented  so  far,  nonlinear 
contributions  are  automatically  introduced  when  we  increase  the  ampli¬ 
tude  of  the  fundamental  mode  (see  Figure  4.5).  This  is  especially  appa¬ 
rent  in  the  second-harmonic  component.  The  amplitude  of  the  second  har¬ 
monics  starts  from  0  at  x  =  1  in  Murdock's  calculation,  while  our  solu¬ 
tion  would  give  non-zero  value  proportional  to  A2.  In  order  to  match 
Murdock's  initial  condition,  we  have  to  construct  unforced  boundary- 
value  solution  for  the  second  harmonic  (2u>  component)  that  will  cancel 
our  forced  second-harmonic  component.  The  free  solution  is  formally 
given  by 


IX  exp  li  k2,n  x  •  21  "t> 

where  and  k0  „  are  the  eigenfunctions  and  the  eigenvalues  of  Orr- 
t,n  Z,n 

Sommerfeld  equation  corresponding  to  the  frequency  2w.  The  expansion 
coefficients  an  are  to  be  determined  from  the  boundary  condition  at 
x  =  1: 


£an  ♦<»,„  .  . 

n 


where  $  (y)  is  the  forced  second  harmonic  at  x  =  1.  Using  the  ortho¬ 

gonality  condition  derived  in  the  appendix,  we  determine 


a_  = 


f 


(2)*  .(2) 
E,n  B  ^E,n 


dy 


(2)* 

where  '  are  the  adjoint  eigensolution  and  B  is  the  operator  defined 
l  ,  n 

by 


B=4k,  (D2-k^  n)  +  iR  [U(D2-3k2  n)  +  4u>k9  -  U‘ '  ] 

2,n  2,n  2,n  2,n 

In  our  calculation,  we  include  only  the  first  mode  since  we  believe  it 
is  the  least  damped  solution.  We  calculated  also  the  correction  to  the 
second-harmonic  eigenvalues  due  to  the  mean-flow  variation. 


X 

k2r 

k2i 

k21  r 

k2 1  i 

2kor 

1. 

65.2 

.28 

m 

1.2 

65.0 

.82 

mfwm 

1.4 

64.8 

1.70 

.41 

-.185 

72.2 

1.6 

64.6 

2.84 

.43 

-.101 

71.8 

1.8 

64.1 

4.17 

.46 

-.004 

71.6 

2.0 

63.4 

5.64 

.49 

+  .135 

71.4 

2.2 

62.3 

7.15 

.46 

.37 

71.3 

2.4 

60.7 

8.39 

.18 

.61 

71.2 

2.6 

58.8 

8.96 

-.23 

.48 

71.1 

2.8 

57.3 

8.96 

-.30 

.19 

71.0 

3.0 

56.2 

8.80 

-.22 

.04 

71.0 

3.2 

55.3 

8.63 

-.15 

-.23 

70.9 

3.4 

54.6 

8.50 

-.09 

-.47 

70.8 

Table  4.6  Second-Harmonic  Eigenvalues 
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We  present  the  computed  eigenvalues  in  Table  4.6.  ^*<2r’*<2i  ^  1S 

eigenvalue  for  2w  frequency  in  local  parallel  flow,  and  ^lc21r,k21i  ^  1S 
the  correction  due  to  the  mean  flow  variation.  The  spatial  amplifica¬ 
tion  rate  of  the  second  harmonics  is  given  by  -^2i+*(21r’  an<*  near  x  =  * 
the  second-harmonic  eigenmode  is  amplified  because  of  the  mean-flow  va¬ 
riation,  while  the  local-parallel-flow  result  is  damped  second  harmo¬ 
nics.  The  last  column  is  twice  the  wave  number  of  the  fundamental  (fre¬ 
quency  w)  eigenmode.  As  we  expect,  k2r  is  close  to  but  not  identical  to 
2kor- 


We  denote  this  free-mode  second  harmonic  component  by 


Re  | A2 ( x )  <j>[2)(y.X)  exp  [ie2(x)  -  2iwt) 

Then,  for  some  distance  downstream  of  x  =  1,  before  it  is  damped  out, 
the  amplitude  of  the  free-mode  second  harmonics  will  be  of  the  same 
order  of  magnitude  as  the  forced-mode  second  harmonics,  and  its  beating 
with  the  basic  fundamental  component  will  produce  an  additional  contri¬ 
bution  to  the  fundamental  component  which  will  be  of  the  form: 

Re  |  A(X)A2(X)4>41)  (y,X)  exp  [i(VV  "  iwt 

The  differential  equation  for  <t>^  is  the  following  nonhomogeneous  Orr- 
Sommerfeld  equation: 


[3y  -Ck2"iro)2]2<j>41)'iR  |t(k2-k0)U(y)-W][3^-(k2-k0)2]4> 


-  (k2-k0)  U  (y)  *4 


(1) 


=  -  iR 


Vty1  (3y-fo>  *<U  +  Me2’  *y'’ 


-  V’11  <’’-kl>  W”  (^-k2>4y) 


(1) 

4 


with  boundary  conditions 


****** 


I 
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♦Jn(0)  =  *^(0)  -  0  and  *  0  as  y+0  . 

Our  results,  including  the  contributions  from  the  free-mode  second  har¬ 
monics,  are  compared  with  Murdock's  computation  in  Figure  4.6.  Since 
the  wave  numbers  for  the  forced  and  free-mode  second  harmonics  are  not 
equal,  the  amplitude  of  the  second  harmonics  is  modulated  with  x  varia¬ 
tion.  Since  we  did  not  include  all  modes  of  the  eigenfunction  for  the 
second  harmonics,  we  do  not  cancel  the  forced  second  harmonics  perfectly 
at  x  =  1,  but  the  second  harmonic  modulation  is  very  well  reproduced 
except  for  a  slight  phase  shift.  Because  of  the  above-mentioned  differ¬ 
ence  in  the  wave  numbers,  the  computed  amplitude  of  the  fundamental  com¬ 
ponent  also  exhibits  small  amplitude  modulation  until  the  free-mode  sec¬ 
ond  harmonics  is  damped  out.  Except  for  this  modulation,  our  results 
for  the  fundamental  component  agree  fairly  well  with  Murdock's  computa¬ 
tions. 

The  longitudinal  component  of  the  perturbation  velocity  at 

nm  =  y-^/ll/2vx  =  0.2  at  a  fixed  t  is  shown  in  Figure  4.7  together  with 

Murdock's  result.  The  agreement  between  the  two  computations  is  good, 
and  especially  the  phase  relationship  between  the  fundamental  component 
and  the  second-harmonics  is  very  well  reproduced  by  the  present  theory. 

From  these  results,  we  conclude  that,  in  the  range  of  x  covered  in 

Murdock's  computation,  the  present  nonlinear  stability  theory,  which 
includes  the  fundamental  and  second-harmonic  components  with  the  correc¬ 
tions  for  the  slow  variation  of  the  undisturbed  laminar  boundary  layer 
flow  with  x,  provide  an  adequate  description  of  the  evolution  of  small 
but  finite  amplitude  perturbations  introduced  into  the  boundary  layer, 
compared  with  direct  numerical  solution  of  the  Navier-Stokes  equa¬ 
tions.  Even  though  the  computation  of  eigensol utions  at  several  sta¬ 

tions  along  the  developing  boundary  layer  Is  quite  involved,  once  the 
eigenvalues  and  eigenfunctions  are  computed,  these  can  be  used  for  any 
amplitude  within  the  limitation  of  the  theory,  and  provide  a  valuable 
"tool"  for  a  parametric  study  in  the  transition  problem. 


9 
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Figure 
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As  an  example,  we  extended  our  computation  to  x  =  3.4  beyond  x  =  2.2, 
and  the  results  are  shown  In  Figure  4.8.  In  this  example,  we  have 
A0  =  18.1  and  the  Instability  burst  occurs  at  x  =  2.73.  It  would  have 
been  very  Interesting  and  valuable  if  Murdock  had  extended  his  computa¬ 
tion  into  this  range  and  obtained  the  flow  field  near  the  Instability 
burst  from  the  numerical  solution  of  the  Navler-Stokes  equations. 


» 


» 


> 


5.  CONCLUSIONS 


We  have  formulated  a  weakly  nonlinear  stability  theory  with  the  assump¬ 
tion  that  the  amplification  rate  is  also  weak.  The  theory  takes  into 
account  the  slow  variation  of  the  undisturbed  mean  flow  with  the  stream- 
wise  distance  and  includes  the  fundamental  and  second-harmonic  compon¬ 
ents.  A  computation  was  carried  out  for  one  frequency  and  the  results 
compare  favorably  with  the  numerical  solution  of  the  Navier-Stokes  equa¬ 
tion  by  Murdock. 

Theory  predicts  that  the  nonlinearity  has  a  stabilizing  effect  at  low 
Reynolds  numbers,  but  becomes  destabilizing  beyond  a  threshold  Reynolds 
number  and  leads  to  a  break-down  of  the  solution.  The  location  of  the 
singularity  depends  on  the  Initial  amplitude  of  the  disturbance,  and 
this  location  may  be  used  as  a  criterion  for  the  laminar- to- turbulence 
transition  replacing  the  commonly  used  en  criterion  which  is  independent 
of  the  initial  amplitude.  The  validity  of  this  new  criterion  should  be 
established  by  carrying  out  the  computation  for  other  frequencies  and 
the  prediction  compared  with  the  transition  location  measured  on  a  flat 
plate  as  the  free-stream  turbulence  level  is  changed. 


The  present  theory  may  readily  be  extended  to  non-Blasius  two-dimen¬ 
sional  boundary  layers. 
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APPENDIX 


Expansion  of  Arbitrary  Function  in  Series  of  Spatial  Mode 
Eigenfunctions  of  Orr-Sommerfeld  Equation 

The  Orr-Sommerfeld  equation  in  the  linear  stability  theory  of  parallel 
flow  is 


(02-k2)2  <j>  -  i  R  C(kU-uj) (02-k2) 4>-kU" 4>]  =  0  (1) 


with  boundary  conditions 

4>  =  D<j>  =  0  aty  =  0 
0  as  y  +  “ 


when  R  and  are  real  constants,  the  problem  is  called  a  spatial -mode 
solution  and  the  eigenvalue  k  may  be  a  complex  number.  Presumably  an 
Infinite  number  of  eigenvalues  (discrete  and  possibly  continuous)  and 
eigensol utions  exist,  and  the  problem  under  consideration  in  this  note 
is  how  one  can  expand  an  arbitrary  function  of  y  in  series  (and  integral 
if  the  eigenvalue  is  continuum)  of  eigenfunctions.  This  problem  was 
considered  and  the  formalism  of  the  expansion  was  given  by  Schensted 
(1961)  for  the  temporal  mode,  in  which  R,  k  are  real  and  is  a  complex 
eigenvalue.  Since  the  coefficients  of  equation  (1)  are  polynomials  of 
the  eigenfunction  k,  the  expansion  formalism  for  the  spatial  mode  is  not 
as  simple  as  in  the  temporal  mode  and,  to  our  knowledge,  it  has  not  been 
given  in  the  literature. 
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Langer  (1923)  states  that  the  systan  of  the  form 

u(n)(x)  +  P1(x,p)u(n_1)(x)  +  ...  +  Pn(x,p)u(x)  =  0 

3* 

with  the  coefficient  P  (x,p)  a  polynomial  of  degree  n  in  p  may  be 
written  in  the  form 

uj(x)  =  E  |pjJt(x)p  +  qjJt(x)i  U4(x) 


E  jo^u^a)  +  |  =  0  ( j=l,2 . n) 

The  Orr-Sommerfeld  equation  falls  into  this  class  of  system. 

Birkhoff  and  Langer  (1923)  developed  the  theory  of  a  system  of  n  ordin¬ 
ary  linear  differential  equations  of  the  first  order  containing  a  param¬ 
eter  and  subject  to  the  homogeneous  boundary  conditions  and,  in  particu¬ 
lar,  discussed  the  formal  development  of  a  vector  of  arbitrary  functions 
into  a  series  of  the  eigensol utions. 

We  are  indebted  to  Professor  D.  Cohen  of  Caltech  (also  a  Consultant  at 
Dynamics  Technology)  for  directing  our  attention  to  the  Langer  and 
Bi rkhoff-Langer  papers  and  for  his  suggestions  in  this  particular  part 
of  our  problem. 


uU)(a) 


u(^(b)  >  =  0 


( j- 1 » 2  *  •  •  • » n ) 
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In  order  to  develop  the  expansion  formalism,  we  consider  the  following 
system  of  first-order  equations  which  is  equivalent  to  the  fourth-order 
equation  given  in  equation  (1).  Let 


z 2  =  4> *  -  k<|>  =  z|  -  kz^ 

z3  =  '  -  k2<i>  =  z'2  +  kz2 

z^  =  <f>' 1 '  -  k<i>'  ‘  -  k24>'  +  k3i)>  =  Z3  -  kz-3 


Then  equation  (1)  becomes 

z^  +  kz4  =  i  R  [ ( kU-oo] ( ^" -k24>]  -  kU  "  <j>] 


Therefore,  equation  (1)  is  equivalent  to  the  system: 
z!  =  kz,  +  z. 


'I 


’I  '  2 

kz~  +  z- 


z3  *  kz3  +  z4 

z^  =  -  iRkU''z^  +  iR(kU-w)z3  -  kz^ 


(2) 


(3) 


The  crucial  feature  of  the  above  system  is  that  the  equations  are  linear 
in  k. 


In  matrix  notation  the  system  (3)  is  written  as 

z’  =  (kA+B)z  (4) 


[*1 ,z2’z3,z4^ 


where 
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A  = 


1 

0 

0 

-  i RU  *  1 


0 

-1 

0 

0 


0 

0 

1 

iRU 


0 

0 

0 

-1 


(4-a) 


B  = 


0  1  0 

0  0  1 

0  0  0 

0  0  - i Rw 

The  boundary  conditions  are 

z,  =  z„  =  0  at  y  =  0  | 


'1  2 
z  +  0 


0 

0 

1 

0 


(4-b) 


as  y  +  »  ( 

We  introduce  the  adjoint  solution  x  =  [x^.x^.x^x^]  that  satisfies 
x'  =  -  x  [kA+B] 

x3  =  x4  =  0  at  y  =  0  ) 

x  +  0  as  y  +  00  ) 


(4-c) 


(5-a) 


(5-b) 


The  reason  for  this  boundary  condition  at  y  =  0  becomes  clear  in  the 
subsequent  development. 


Now,  suppose  that  z^  corresponds  to  an  eigenvalue  km  and  x_n  corresponds 
to  another  eigenvalue  kn;  namely 


zL  =  ( k  A+B)  z 
~m  m  ~m 


'  -  Jn  ,knAtB) 
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By  cross-multiplying  and  adding  the  products,  we  obtain 


( x  z  )'  =  (k-k  )  xn  A  z 
~n  -sm  m  n  ~n  ~m 


By  integrating  with  respect  to  y  from  0  to  00 ,  we  get 


“vk 


.>J\ 

»  A 


A  L  dy  =  x„  z 
~m  [_~n  ~m  J 


xnlzml  +  xn2zm2  +  xn3zm3  +  xn4zm4 


The  integrated  terms  vanish  by  the  boundary  condition  imposed  on  x. 
Therefore,  we  obtain  the  ortho-normality  condition 


x  A  z  dy  =  6 
~n  •'in  3  nm 


(6) 


Jo 


Eigenfunction  Expansion; 


If  an  arbitrary  function  z  is  expanded  in  series  of  z^ 
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Hence 


/•» 


am  = 


x  A  z  dy 
~tn  ~  J 


(8) 


Transformation  to  Single-Equation  Formalism: 

The  relation  between  the  solution  of  equation  (1),  <P,  and  z  is  given  in 
equation  (2).  Written  for  each  component,  equation  (5)  is 


xi  =  kxi  +  ^u"  x4 

x2  =  "  X1  +  kx2 
X3  =  '  X2  “  kx3  “  -’R(kU-u)x4 
-  -  X3  +  kx^ 

By  eliminating  x1,x2,x3,  we  get 

(D2-k2)Z  x4  -  D2  [fR(kU-w)x4]  -  iR  [k2(kU-«]+kU' ' ]  x4  =  0 


wi  th 


at  y  =  0 
at  y  =  ” 


This  is  the  adjoint  equation  directly  obtainable  from  equation  (1).  If 
we  denote  the  adjoint  solution  by  $*,  we  get  the  following  transforma¬ 
tion: 
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x 

x 

X 

X 


s*  -D 3 4>*  +  kD2<J>*  +  D  [k2  +  iR(kU-u))]«j>* 
2  =  D2<Ji*  -  [k2+iR(kU-w)]4>* 

=  _  d«*  +  k$* 


-  [k3  +  iRk(kll-u))]^* 


If  f(y)  is  the  first  component  of  a  function  z  to  be  expanded,  then  the 
other  components  are  given  by  the  relation  in  equation  (2): 


(  Zl  =  f 

1  z2  =  Df  -  kf 

|  y  D2f  -  k2f 

(  z4  =  D3 f  -  kD2f  -  k2Df  +  k3f 


Therefore,  the  formula  given  by  equation  (8)  becomes 


0 


am  =  I  \  (4k  +iRU)  D2f  -  [4k3+iR(3k2U-2k  u+U'  ’ )  ]  f  }  <t>*  dy 
mllm  m  m  m  lmJ 


The  normalization  relation  is 


<4V™1  °2** 


[4k>iR(3kjlJ-2k  ' )]» 
m  mm  m 


dy  «  1 
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